#!/usr/bin/env python
# coding: utf-8

# In[1]:


import pandas as pd
import numpy as np
import os
import warnings
from pandas.tseries.offsets import Hour,Minute,Day
warnings.filterwarnings('ignore')


# ## 读取数据集

# In[2]:


import json
with open('flight_dict.json','r') as f:
    flight_dict = json.load(f)
product = pd.read_csv('data_fam_products.csv')
sche = pd.read_excel('time_info.xlsx')
fleet = pd.read_csv('data_fam_fleet.csv')


# ## 数据预处理

# ### 处理产品信息

# In[3]:


product.head()


# In[4]:


product['TRD'] = product.loc[:,product.columns[product.columns.str.contains('RD')]].sum(axis=1)


# In[5]:


product = product.drop(product.columns[product.columns.str.match(r'RD\d')],axis=1)


# In[6]:


product.head()


# In[7]:


product = product[product.TRD != 0]
product.describe()


# In[8]:


def f(x):
    if x != '.':
        return x[:6]
    else:
        return x


# In[9]:


product.Flight1 = product.Flight1.apply(f)
product.Flight2 = product.Flight2.apply(f)
product.Flight3 = product.Flight3.apply(f)


# In[10]:


all_flight = pd.unique(pd.concat([product.Flight1,product.Flight2,product.Flight3]))
all_flight = all_flight[all_flight != '.']
print(len(all_flight))


# In[11]:


product.head()


# In[12]:


product[product.Flight1.str.contains('AA0040')]


# In[13]:


product.index = range(len(product))


# In[14]:


TRD_ = product.TRD
Fare_ = product.Fare
TRD_

import json
flight_dict = {}
for flight in all_flight:
    flight_dict[flight] = list(product[(product.Flight1.str.contains(flight)) | \
                                  (product.Flight2.str.contains(flight)) | \
                                  (product.Flight3.str.contains(flight))].index)
with open('flight_dict.json','w') as f:
    json.dump(flight_dict,f,indent=2)                
# In[15]:


flight_table = pd.Series(flight_dict)
flight_table = flight_table.sort_index()
flight_table.head()


# ### 处理时刻表信息

# In[16]:


sche.head()


# In[17]:


sche.dtypes


# In[18]:


sche.dep_ucl.min()


# In[19]:


sche.arr_ucl.min()


# In[20]:


sche[sche.arr_ucl ==sche.arr_ucl.max()]


# In[21]:


sche


# In[22]:


sche.flight_time = sche.flight_time.astype(str)


# In[23]:


def cacl_time(x):
    h,m,s = x.split(':')
    return int(h)+int(m)/60


# In[24]:


sche.flight_time = sche.flight_time.apply(cacl_time)


# In[25]:


sche.head()


# In[26]:


sche_g = sche.groupby('flight')


# In[27]:


every_flight_time = sche_g.flight_time.sum()


# In[28]:


every_flight_max_TRD = sche_g.TRD.max()


# In[29]:


every_flight_min_dep_ucl = sche_g.apply(lambda x:x.loc[x.dep_ucl == x.dep_ucl.min(),['origin','dep_ucl']])
every_flight_max_arr_ucl = sche_g.apply(lambda x:x.loc[x.arr_ucl == x.arr_ucl.max(),['destination','arr_ucl']])


# In[30]:


every_flight_max_arr_ucl = every_flight_max_arr_ucl.reset_index().drop(['flight','level_1'],axis=1)
every_flight_max_arr_ucl.head()


# In[31]:


every_flight_min_dep_ucl = every_flight_min_dep_ucl.reset_index().drop(['level_1'],axis=1)
every_flight_min_dep_ucl.head()


# In[32]:


every_flight_max_TRD = every_flight_max_TRD.reset_index().drop(['flight'],axis=1)
every_flight_max_TRD.head()


# In[33]:


every_flight_time = every_flight_time.reset_index().drop(['flight'],axis=1)
every_flight_time.head()


# In[34]:


final_df = pd.concat([every_flight_max_arr_ucl,every_flight_min_dep_ucl,every_flight_max_TRD,every_flight_time],axis=1)


# In[35]:


final_df.head()


# In[36]:


final_df.dep_ucl.min()


# In[37]:


final_df.arr_ucl.max()


# In[38]:


final_df[['arr_ucl','dep_ucl']] -= final_df.dep_ucl.min()


# In[39]:


final_df.head()


# In[40]:


final_df = final_df.reindex(columns=['flight','origin','destination','dep_ucl','arr_ucl',
                         'flight_time','TRD'])


# In[41]:


final_df.head()


# In[42]:


final_df[[f'机型{i}' for i in range(9)]] = fleet.每小时飞行成本.values * final_df.flight_time.values[:,np.newaxis]
cost_df = final_df.loc[:,['flight']+list(final_df.columns[-9:])]
cost_df = cost_df.set_index('flight')
cost_df.head()


# In[43]:


fleet.index = final_df.columns[-9:]
fleet.head()


# In[44]:


fleet


# In[45]:


final_df.head()


# In[46]:


final_df.arr_ucl += Minute(40)


# In[47]:


final_df[(final_df.dep_ucl >= Day(1)) & (final_df.arr_ucl >= Day(1))]


# In[48]:


inter_flying = final_df[(final_df.dep_ucl < Day(1)) & (final_df.arr_ucl >= Day(1))].flight
inter_flying


# In[49]:


final_df.arr_ucl[final_df.arr_ucl >= Day(1)] -= Day(1)
final_df.dep_ucl[final_df.dep_ucl >= Day(1)] -= Day(1)


# In[50]:


final_df[(final_df.dep_ucl <= Day(1)) & (final_df.arr_ucl >= Day(1))]


# In[51]:


final_df.loc[:,['dep_ucl','arr_ucl']].min()


# In[52]:


final_df.arr_ucl = (final_df.arr_ucl.dt.seconds / 60).astype(int)
final_df.dep_ucl = (final_df.dep_ucl.dt.seconds / 60).astype(int)


# In[53]:


final_df.head()


# In[54]:


time_arr = np.sort(pd.unique(pd.concat([final_df.dep_ucl,final_df.arr_ucl])))


# In[55]:


time_dict = dict(zip(time_arr,range(len(time_arr))))


# In[56]:


final_df.dep_ucl = final_df.dep_ucl.map(time_dict)
final_df.arr_ucl = final_df.arr_ucl.map(time_dict)


# In[57]:


len(time_arr)


# In[58]:


final_df.head()


# 生成每个地点的离开时刻集合

# In[59]:


origin = final_df.groupby('origin')


# In[60]:


from collections import defaultdict
group_dep = {}
for label,group in origin:
    group = group.sort_values(by='dep_ucl')
    group_dep[label] = defaultdict(list)
    for i in range(len(group)):
        dep_t = group.iloc[i].dep_ucl
        flight = group.iloc[i].flight
        group_dep[label][dep_t] += [flight]


# 生成每个地点的到达时刻集合

# In[61]:


des = final_df.groupby('destination')
group_arr = {}
for label,group in des:
    group = group.sort_values(by='arr_ucl')
    group_arr[label] = defaultdict(list)
    for i in range(len(group)):
        arr_t = group.iloc[i].arr_ucl
        flight = group.iloc[i].flight
        group_arr[label][arr_t] += [flight]


# In[62]:


group_arr


# In[63]:


group_dep


# In[64]:


pd.concat([pd.Series(group_arr['BER']),pd.Series(group_dep['BER'])],axis=1)


# In[65]:


L_dict = {}
Node_dict = {}
origin_label = set(origin.groups.keys())
des_label = set(des.groups.keys())
all_label = origin_label | des_label
for label in all_label:
    dep_time = pd.Series(group_arr[label])
    ori_time = pd.Series(group_dep[label])
    both_time = pd.concat([dep_time,ori_time],axis=1)
    both_time.columns = ['in','out']
    L_dict[label] = both_time
    Node_dict[label] = len(both_time)+1 


# ## 开始计算

# In[66]:


import gurobipy


# In[67]:


Model = gurobipy.Model('The plane distribution')


# ## 添加变量

# ### 添加x变量

# In[68]:


x = Model.addVars(cost_df.shape[0],cost_df.shape[1],lb=0,
              vtype=gurobipy.GRB.BINARY,name='x')
x = pd.DataFrame(np.array(x.select('*','*')).reshape((-1,9)),
                 index=cost_df.index,columns=cost_df.columns)
x.head()


# ### 添加y变量

# In[69]:


y_array = np.zeros(9)
for label,value in Node_dict.items():
    y_ = Model.addVars(value,9,lb=0,vtype=gurobipy.GRB.INTEGER,name=f'{label}_y')
    y_array = np.vstack([y_array,np.array(y_.select('*','*')).reshape((-1,9))])
y_array = y_array[1:]


# In[70]:


Node_num = np.array(list(Node_dict.values()))
cumsum_Node = Node_num.cumsum()
y_result_num = []
for i,label in enumerate(Node_dict.keys()):
    ls = list(range(Node_num[i]))
    for j in ls:
        y_result_num.append((label,j))


# In[71]:


y = pd.DataFrame(y_array,index=pd.MultiIndex.from_tuples(y_result_num),columns=x.columns)


# In[72]:


y.head()


# ### 添加z变量

# In[73]:


product


# In[74]:


z = Model.addVars(len(product),1,vtype=gurobipy.GRB.CONTINUOUS,name = 'z').select('*','*')
z = pd.Series(z)
z


# In[75]:


Model.update()


# ## 设定目标

# In[76]:


flight_num= fleet.飞机数量
fleet_values = fleet.座位数
fleet_values.index = x.columns
flight_num.index = x.columns


# In[77]:


s = (z*Fare_).sum()-(x*cost_df).sum().sum()
Model.setObjective(s,gurobipy.GRB.MAXIMIZE)


# ## 添加约束

# ### 约束1：每个航班飞机数量不超过1

# In[78]:


Model.addConstrs((gurobipy.quicksum(x.iloc[i,:]) <= 1 for i in range(len(x))),
                    name='航线飞机数量约束');


# ### 约束2：每个机场的每个时间节点流入等于流出

# In[79]:


contra_in_out = {}
for label in Node_dict.keys():
    time_Node = L_dict[label]
    contra_in_out[label] = []
    for i in range(len(time_Node)):
        in_,out_ = time_Node.iloc[i,:]
        x_in = np.zeros(9)
        x_out = np.zeros(9)
        if isinstance(in_,list):
            x_in = x.loc[in_,:].sum(axis=0)
        if isinstance(out_,list):
            x_out = x.loc[out_,:].sum(axis=0)
        y_in = y.loc[label,i]
        y_out = y.loc[label,i+1]
        contra_in_out[label].append(x_in + y_in - x_out - y_out)
        Model.addConstrs((x_in[j] + y_in[j] - x_out[j] - y_out[j] == 0 for j in range(9)),name=f'constr_{label}_{i}')


# ### 约束3：每个机场的回旋流约束

# In[80]:


for label in Node_dict.keys():
    Model.addConstrs((y.loc[label].iloc[0,i] == y.loc[label].iloc[-1,i] for i in range(9)),name=f'{label}_回转流');


# ### 约束4：初始点设定约束

# In[81]:


init_values = x.loc[inter_flying].sum(axis=0) + y.swaplevel().loc[0].sum()


# In[82]:


Model.addConstrs((init_values[k] <= flight_num[k] for k in range(9)),name='飞机数量约束');


# ### 约束5：飞机座位约束

# In[83]:


flight_table


# In[84]:


seat_values = fleet_values
for flight in flight_table.index:
    Model.addConstr((seat_values*x.loc[flight]).sum() - (z.loc[flight_table[flight]]).sum() >= 0,name='座位约束')


# ### 约束6：实际供应尽可能不大于需求

# In[85]:


Model.addConstrs((z[i] <= TRD_[i] for i in range(len(TRD_))),name='供给约束');


# ## 模型参数设定和训练

# In[ ]:


from gurobipy import GRB,GurobiError
try:
    Model.Params.TimeLimit = 3600
    Model.setParam(GRB.Param.PoolSolutions, 1024)
    Model.setParam(GRB.Param.PoolGap, 0.3)
    Model.setParam(GRB.Param.PoolSearchMode, 2)
    Model.update()
    Model.optimize()
except GurobiError:
        print('Error reported')


# In[ ]:


status = Model.Status
status


# In[ ]:


GRB.OPTIMAL


# In[ ]:


Model.getObjective().getValue()


# In[ ]:


import joblib
joblib.dump(Model,'current_model.pkl')


# In[ ]:


Model.write('current_model.lp')
Model.write('current_model.mps')
Model.write('current_model.sol')


# In[ ]:


nSolution = Model.SolCount
print('number of solution stored:' + str(nSolution))

# Print objective values of solutions
for e in range (nSolution):
    Model.setParam(GRB.Param.SolutionNumber, e)
    print ('%g' % Model.PoolObjVal, end=' ')
    if (e + 1) % 5 == 0:
        print('')


# In[ ]:


for k in range(9):
    print(init_values[k].getValue())


# In[ ]:


result = Model.getVars()
x_result = pd.DataFrame(np.array(result[:6651]).reshape((-1,9)),columns=x.columns,index=x.index)
x_result = x_result.applymap(lambda x1:x1.x).astype(int)
x_result


# In[ ]:


x_result.sum()


# In[ ]:


len(z)


# In[ ]:


y_result = result[6651:-20053]
y_result =pd.DataFrame(np.array(y_result).reshape((-1,9)),columns=x.columns)
y_result = y_result.applymap(lambda x1:x1.x).astype(int)
y_result


# In[ ]:


z = z.apply(lambda x1:x1.x)


# In[ ]:


z


# In[ ]:


y_result = pd.DataFrame(y_result.values,index=pd.MultiIndex.from_tuples(y_result_num),columns=y_result.columns)


# In[ ]:


y_result


# In[ ]:


y_result[y_result.机型3 == 2]


# In[ ]:


product['z'] = z
y_result.to_excel('y_result.xlsx')
x_result.to_excel('x_result.xlsx')
product.to_excel('z_result.xlsx')

